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Abstract 

A numerical  finite  difference  method 
is  developed  to  simulate  the  viscous  flow 
over  blunt/sharp  bodies  at  incidence. 
Herein,  a two-layer  model  Is  suggested. 
The  inner  region  consists  of  the  three- 
dimensional  boundary  layer  and  boundary 
region . Laminar  and  turbulent  flows  are 
considered.  The  governing  system  applies 
in  boundary  regions  and  for  problems  with 
cross  flow  reversal.  The  equations  are 
integrated  by  a predictor-corrector 
scheme.  For  the  turbulent  boundary  layer 
analysis,  both  a mixing  length  model  and 
a two-equation  kinetic  energy-dissipation 
system  is  considered  for  Reynolds  stress 
closure.  The  outer  region  is  inviscid. 
The  Euler  equations  are  integrated  with 
McCormack's  two  level  explicit  scheme. 

For  the  matching  of  the  two  regions, 
three  dimensional  viscous  displacement 
and  entropy  layer  swallowing  are  consid- 
ered. Numerical  solutions  are  compared 
with  experimental  data  and  indicate  that 
the  present  formulation  can  give  an  ac- 
curate prediction  of  aerodynamic  loads, 
skin  friction  and  heat  transfer  rates  on 
sphere-cone-cylinder-flare  shape  bodies 
at  angle  of  attack.  The  calculations  are 
suitable  for  (i)  supersonic  or  hypersonic 
freestreams,  (11)  large  Reynolds  number 
and  (iii)  flows  without  streamwise  flow 
separation;  however,  secondary  flow  re- 
versal is  allowed. 


F = U/Uc 


Nomenclature 

u,  V, 

x,  y. 

a , b,  k 

constants  defined 
eq.  (24) 

in 

c 

constants  defined 
eq.  (14) 

in 

II 

A = 26  Jpw/s> 

d = <5-  & 

A 

E,  F,  G,  I 

Matrices  defined 

in  eq.  (1) 

6 

h constant  defined  in  eq - (13) 

or  local  enthalpy 

hv,  hA  metric  coefficients  defined 

' in  eq.  (lla) 

H total  enthalpy 

K turbulence  kinetic  energy 

L characteristic  length 

M Mach  number 

p pressure 

Pr  Prandtl  number 

- („2  ♦ »V/2 

a surface  heat  transfer 

V 

<^cw  = Tr/To  " Tw/ToI 

Re0  = 9/ S'*. 

J surface  length 

t =J  A f J ^ ^ 

T static  temperature 

u,  v,  w velocity  components 

x,  y,  z coordinate  frame 

ck.  angle  of  attack 


variable  defined  in  eq . (6) 


three  dimensional  viscous 
displacement  thickness 

boundary  layer  thickness 


-ff  ('-tity, 


♦This  research  is  jointly  sponsored  by  the  Air  Force  office  of  Scientific  Research 
(under  Grant  No.  AFOSR  7^-2635)  and  the  Avco  Independent  Research  and  Development 
program . 


~£‘  {'-s’/*  <*)  •>  a 

6 local  body  slope 

A = 4i  A*  Aj 

Al  variable  defined  in  eq.  (24) 
turbulent  Eddy  viscosity 
molecular  viscosity 


coordinate  frame 


shear  stress  at  the  wall 
density 

turbulence  dissipation  function 
constant  defined  in  eq . (17) 


subscripts 


refe-rs  to  shock 


b,  w refers  to  surface  condition 


recently  turbulent  flow.'-'  However, 
these  efforts  are  still  largely  confined 
to  either  two  dimensional  or  axisym- 
metric  geometries  and  even  then  are  quite 
time  consuming.  A second  approach  uses 
a simplified  Navier-Stokes  system.  Its 
salient  feature  is  a single-layer  mod- 
el( 10-13)  that  also  precludes  the  neces- 
sity for  matching  of  the  inner  and  outer- 
regions.  Although  this  approach  has  been 
successful  in  a number  of  problems,  the 
computation  times  can  become  excessive. 
One  inherent  difficulty  in  the  single 
layer  model  is  the  specification  of  one 
coordinate  system  to  account  for  the  two 
different  length  scales  in  the  vl  ecus 
and  lnviscid  regions.  For  example  the 
laminar  boundary  layer  thickness  may  in- 
crease as  fyT  (x  - streamwise  distance), 
while  the  outer  bow  shock  grows  linearly. 
It  is  difficult  for  a single  coordinate 
transformation  to  take  into  considera- 
tion the  different  growth  rates  associ- 
ated with  these  two  distinct  regions. 
Therefore,  the  two-layer  model  proposed 
here  appears  to  be  a practical  and  eco- 
nomical alternative  for  weak  interaction 
three-dimensional  flow  calcul -ticn.'  . 


e refers  to  boundary  layer  edge 
conditions 

indicates  b (■ 
respectively 

00  refers  to  free  stream  conditions 


I.  Introduction 

During  the  last  decade,  numerical 
computations  for  both  boundary  lay- 
ersl-'-"3)  and  invlscid  flows'1*-?)  have 
been  advanced  considerably  in  capability 
speed  and  accuracy.  However,  the  com- 
putational procedures  for  these  two  flow 
regions  have  been  primarily  developed 
independently  so  that  the  effects  of 
vlscous-inviscld  interactions  have  been 
largely  neglected. 

In  many  cases  an  lnviscid  calcula- 
tion may  become  inaccurate  unless  proper 
consideration  is  given  to  the  surface 
viscous  interaction.  On  the  other  hand, 
accurate  boundary  layer  edge  properties 
are  sometimes  difficult  to  obtain  unless 
a meaningful  lnviscid  computational 
program  is  available.  Therefore,  numer- 
ical calculations  of  the  inner  (vis- 
cous) and  outer  (lnviscid)  regions 
should  be  properly  interrelated  by  suit- 
able matching  conditions  . While  this 
has  been  done  for  two  dimensional  flows, 
three  dimensional  matched  solutions  for 
general  geometries  have  not  been  pre- 
viously obtained. 

Direct  integration  of  the  complete 
Navier-Stokes  equations  avoids  the 
matching  process.  This  has  been  car- 
ried out  successfully  byqa  number  of 
researchers  for  laminar' °>  and  more 


The  formulation,  governing  equations 
and  numerical  methods  for  the  inner 
(viscous)  and  outer  (lnviscid)  regions 
will  be  briefly  discussed  in  sections  II 
and  III,  Typical  numerical  solutions  are 
compared  with  existing  experimental  data. 
The  Inner  and  outer  solutions  are  then 
coupled  by  Including  the  pertinent  ef- 
fects of  viscous  displacement  and  en- 
tropy layer  swallowing  effects.  This 
procedure  is  described  In  section  IV 
where  some  relevant  examples  are  depict- 
ed . 


II,  Outer  lnviscid  Flow 

The  three  dimensional  lnviscid  flc-w 
in  the  shock  layer  of  a blunt  non-cir- 
cular body  at  angle  of  attack  is  consid- 
ered. A time  dependent  finite-differ- 
ence technique  is  applied  for  the  sub- 
sonic nose  region  and  streamwise  march- 
ing procedure  for  the  supersonic  after- 
body. A two-level  explicit  scheme  is 
used  for  the  calculation  of  the  internal 
points  and  a modified  method  cf  charac- 
teristics procedure  is  applied  at  bound- 
ary and  shock  points. 


II. 1 Supersonic  Region 


The  lnviscid  formulation  is  patter- 
ned after  that  of  Moretti ; ' *** » how- 

ever, the  governing  equations  are  writ- 
ten in  divergence  form.  The  Euler  equa- 
tions in  cylindrical  coordinate  are: 


£ y +■  Gq  +■  Z - O 


and  more 


2 


where 


pit*  + p/m£  s' 

£ ~ Pti  v 

Pa 

P W r'1 

Pawr~‘ 

Pvw  r -i 

(pw*+  7/MJrjr-1: 

f - e 7^  ^ +pVuv 

V/M^+PV* 

f>vw 

P vr 

I ~ Puvr'^ 

P(v2  - w2)  r '*  + A,  E + a2& 
P vw 


A-A<^Z/~  7Z 

+ k J[  ( 7*)  + 7r  (^'k)f<Y 

Ufa  +-  nr+p7*>^K'tiiJ' 

McCormack's  :hem<  ' a]  i to  ' nfc  - 

grate  eq.  ' for  (i  r } cy^r- 

is  evaluated  ty  a three-point  end  dif- 
ference formula.  It  should  tv  noted  that 
iteration  is  unnecessary  fcr  tiv  form- 
ulation cn  boundary  points. 

On  the  body  surface,  it  is  require! 

that  — r 

f v rt  = o 

(5) 

The  wall  pressure  is  obtained  from  the 
compatibility  relation  along  the  left 
running  characteristic,  i.e, 

P?=  -a  P?  h A\UV  v -( VPr^)  Ul\ 

+ U Vl  + aW‘  (rfy  + ? 

//  r r 


A = [(H),  -&)*]/[ r,  - ri.] 

7 = (r-ny  ft  - rt) 

Here  U,V,  Wj  P and^3 

are  nondimensionalized  with  respect  to 
free  stream  flow  properties . The  temp- 
erature distribution  is  evaluated  from 
the  following  isoenergtic  conditon, 

T--v/p--\+  »:( *■')(  uvwj) 

An  explicit  numerical  method  is 
employed.  The  two  level  integration 
scheme  suggested  by  McCormack'1®'  has 
been  adopted  for  the  interior  points. 

The  marching  A 2 is  limited  by  the 
Courant-Frederick-Levy  conditon. 

The  outer  bow  shock,  r = r ( 2 , <P  ) 
is  treated  as  discontinuity®  Flow  prop- 
erties behind  the  shock  are  evaluated 
from  the  Rankine-Hugoniot  relation  once 
the  normal  velocity  component 
0 • v r is  obtained . An  additional 
condition  which  determines  the  value  or 
equivalently  the  shock  orientation  is 
supplied  by  the  following  characteristic 
compatibility  relation: 

**  rs/ f(Sr>  *V)  (3) 

The  complete  expression  for 
is  derived  in  Reference  (15).  Eq . (3)  is 
valid  along  the  right  running  character- 
istic surface  defined  by 


Since  the  body  entropy  is  known,  the  den- 
sity can  be  computed  from  the  lsentropic 
relation.  Finally  the  velocity  compon- 
ents at  wall  are  calculated  from  eqs . 

(2)  and  (5). 

11. 2 Invlscld  Blunt  Body  Calculation 

( u ) 

Moretti's'  ’ time-dependent  blunt 
body  program  has  been  used  tc  calculate 
the  flew  field  in  the  subsonic  nose  re- 
gion. This  code  is  mere  versatile  than 
the  inverse  method  of  Pef.  (7),  particu- 
larly when  the  free  stream  Mach  numb< i Is 
lew  (Moc=?),  and  the  body  is  oth  r then 
the  sphere-cone  geometry , Recently  this 
numerical  program  has  been  modified  tc 
include  real  gas  and  nonuni  form  fre 
stream  effects.'  ' The  blunt  tody  sol- 
utions supply  the  starting  renditions  for 
the  downstream  supersonic  flew  calcula- 
tions . 

11. 3 Results  fcr  Invlscld  Flow 

In  order  to  check  the  computer  code 


3 


III. 


Formulation 


> 


and  to  demonstrate  Its  capability  several 
sample  calculations  are  presented  here. 

The  surface  pressure  distribution  over 
an  elliptic  cone  (b/A  - 1.79)  at  15° 
angle  of  attack  Is  depicted  In  Fig.  2. 
Comparison  is  made  with  Zakkay  and 
Visich's  data:'-1-')  the  agreement  is  en- 
couraging. It  should  be  noted  that  for 
ot/e>  1.5  ( ct-  = angle  of  attack,  9 = 
cone  half  angle)  the  inviscid  computa- 
tional for  a sharp  cone  becomes  unstable 
and  the  inner  viscous  region  cannot  be 
neglected.  The  aim  of  our  present  work 
is  to  correct  the  deficiency  of  the  in- 
viscid solution  by  matching  with  the 
boundary  layer. 

Fig.  3 shows  the  inviscid  results 
over  a sphere  cone  (9  = 9°)  at  angle  of 
incidence  ( oL  = and  10°).  The  ex- 
pansion-recompresslon  processes  near  the 
junction  of  sphere  and  cone  are  simulated 
quite  well  by  the  finite  difference  sol- 
utions . 

The  last  case  considered  is  the  sur- 
face pressure  distribution  on  a cone- 
cylinder-flare  which  is  shown  in  Fig.  4. 
Also  presented  are  Zakkay  and  Callahan's 
experiment  measurements .v 1°) 

In  summary,  it  has  been  demonstrated 
that  the  numerical  code  for  the  inviscid 
flow  is  reasonably  accurate  and  versatile. 
In  the  following  section,  the  formula- 
tion for  the  inner  region  is  discussed. 


III.  Inner  Layer  (Viscous  Region) 

Many  of  the  existing  theoretical 
studies  of  three  dimensional  boundary 
layer  either  assume  a similarity  approxi- 
matlon(19)  (such  that  the  governing  eq- 
uations become  pseudo  two  dimensional), 
or  resort  to  a small  cross  flow  approxi- 
mation'20, 21 ) or  independence  principle. 
The  nature  of  the  three  dimensional  bound- 
ary layer  equations  has  been  investigated 
by  Der  and  Raetz1.  Recently  Davis'2^’  c4) 
and  his  associates  have  made  a series  of 
studies  on  numerical  methods  for  nonslmi- 
lar  boundary  layers  on  sharp  as  well  as 
blunt  bodies.  Generally  the  Dwyer-Krause 
method(25>  26)  is  suitable  for  boundary 
layer  computations,  but  without  the  exist- 
ence of  boundary  regions.  The  formula- 
tion suggested  by  Lin  and  Rubin'2'~30) 
can  handle  problems  with  boundary  regions. 
One  advantage  of  this  method  is  the 
ability  to  resolve  flows  with  crossflow 
reversal.  Significantly  this  procedure 
also  removes  the  difficulty  concerning 
existence  and  uniqueness  of  the  solutions 
near  the  leeplane . ' 30a ) 

Recently,  Blottner  and  Ellls^ 
have  generalized  the  Dwyer-Krause  scheme 
for  laminar  incompressible  flow,  further- 
more, they  suggest  a system  of  useful 
coordinates  for  blunt  body  calculations. 


It  ha.  b :n  aosarycd  experimentally-'’2’ 1 
and  analytically' 1 ‘ » 2 - > that  s-.condary' 
flow  reversal  does  not  occur  at  the  t:p  of 
a sharp  cone,  or  in  a blunted  nose  region 
when  the  body  is  at  moderate  angle  of  in- 
cidence. This  important  observation  is 
Implicit  In  the  present  theoretical  form- 
ulation. Herein,  the  Blcttnor  and  Eli:, 
procedure (31)  will  be  adopted  for  flow  in 
the  blunt  no  e region,  while  the  predictor- 
corrector  formulation' ^ '"' °)  is  used  for 
the  afterbody  supersonic  flew . This  choice 
is  acceptable  as  a boundary  region  does 
not  appear  n-.-ar  the  nose,  but  cross  flov. 
reversal  is  possible  in  the  downstream 
flow.  These  two  approaches  were  original- 
ly designed  for  laminar  flow,  but  are  ex- 
tended to  turbulent  flow  conditions  in  the 
present  paper. 

III.l.A  Blunt  Body  Region 

A body  orientated  coordinate. system 
suggested  by  Blottner  and  Ellis'-'-1)  Is 
employed  here  (Fig.  5).  One  coordinate 
is  defined  from  the  intersection  cf  the 
body  surface  with  the  plane  containing 
the  x-axis  and  at  an  angle  <p  from  the 
y-axis.  These  lines  define  the  coordinate 
0 = constant,  (here  one  would  choose 
the  y-axis  to  be  the  symmetry  plane).  The 
other  coordinates  are  orthogonal  4 = con- 
stant lines,  (see  Fig.  5). 

The  surface  of  the  body  is  defined  by 
an  expression  cf  the  form 

"i  - X,  (*z,<p)  , X2r  X2  (X,,  4) 

where  x and  x _ represent  x and  r 
respectively.  The2position  vector  r can 
be  written  as 

r = x i + rcos4-|  + rsin^  k 

Vectors  that  are  tangent  to  = constant 

and  4*  = constant  areJ^^J^  and£=(* 

respectively.  For  the  coordinates  to  be 
orthogonal,  the  relation  s»t  = 0 must  be 
satisfied.  This  leads  to 

Ax,\  _A . (j<rb.  on 

Along  4 = constant,  this  relation  is  writ- 
ten in  finite-difference  form  as 

=(*i)«+  AK  + fA  (8) 

For  points  along  the  symmetry  line,  we 
obtain: 

AX*  =A Slfl  + (AXi/a 

1 ' ' (9) 

A S = hfdi  along  <f>  = 0 


■ 


U 


The  metric  coefficients  are  found  from 
the  definition: 

y / 

+ driJ  I d 4 


(10) 


/ty=|Jx*-f-<yr*+(rd4>)z|  / d <}> 


Eqs  „ (10)  are  expressed  in  finite  differ- 
ence form  in  the  numerical  computations . 

If  the  inviscid  flow  on  the  body 
surface  is  given  in  terms  of  the  (x,  y, 
z)  coordinates  (Pig.  5)>  so  that 

UT=  U,  I 4-  Vi  ^ + W,  te. 

then  the  velocity  components  parallel  to 
the  £ and  4>  planes  can  be  found  from 

Ue  = (Ur  ■ 5)/ IS/ 

1*4  — (llr  • t) / 1 1 I 

The  three  dimensional  boundary  lay- 
er equations  in  terms  of  the  curvilinear 
coordinat  s ) can  be  written  as: 

continuity 
4 momentum 

* v*  tlf|-  Hf*K-  * 

-I  327 


■ v n />? 

a-  momentum 


W, 

d 


7 + { f wr9  *i F w ~ w^|  + 

JuAi i = _ _L  Aiit/Fz-®)t/w*- 

d$  I a 4>  ' ' ' 


_&y 


(ID 


energy 


+ Tri ' V r-i ) 

(*1  + ^7)  * Tr  * 7T  ®( F "W  &) 

where 

F=  ^6.  , U/  = ,®= 


A~  j ~ - ^A 

The  metric  coefficients  ht  and  1 4,  which 
were  given  in  eq.  (10)  are  defined  ‘a.: 


dt-  +■  d$  t dy‘ 


(lla) 


Along  the  symmetric  line  these  equations 
can  be  simplified  to: 


continuity 
y momentum 

v„v^|FFs+y- «)  y*u£\^&f7) 

momentum  ' 

Ve  6?  + { F&i  H F&  - ®J 

+ ^Jjliltl  4- -(4-  (62-®)-(-  A We  _ a 
dt  I hit  A 


energy 


i<t> 


(12) 

hfisMu 


MH  ^(^7)  ’ ^ 

where  6=  W/W* 

Near  the  stagnation  point  (^  =0),  Taylor 

series  expansions  have  been  used  for  W , 

Ue,  and  h*  , 


We  = al£ 


Ue  = bl£ 


h = h £ 


(13) 

It  should  be  noted  that  a distinction’  1 ^ 
must  be  made  between  the  cases  when 
a.  ^ 0 and  a.  = 0.  Eq.  (11)  can  be  fur- 
ther simplified  to  the  following  form  at 
the  stangation  point, 

continuity 

(^\+F+iVl6f7^6^r(|r)|+6i7=^ 

£ momentum 
when  Cb  t 0 

14  + F6  -®  + A-  ±(<; 
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when  a-^  = 0, 


+ (Ft  -®)-f  (f-®)=  f ( 

£ momentum 

Vi  ff+-2&  fj-  +-£  {F  -®) 

+ffc(rc-®)l&  $-1 

-F®>-®)=£r(^7) 

energy:  ' 


There  is  no  restriction  on  the  magnitude 
of  w in  eq.  (11-14),  therefore,  the  gover- 
ning system  is  not  limited  to  small  cross- 
flow.  Finally  the  boundary  conditions 
for  eqs . (11-14)  are 

Tj=.0  , F-  W -Vt  - & - O 
® = ®w  or  (g^j  - 0 

7 = V , f-  ® = & = / 

W = Wfc/Ue 

II1.1.B  Afterbody  Supersonic  Region 

Downstream  of  the  blunt  nose,  the 
formulation  described  in  Ref.  30  is  ad- 
opted here.  The  conventional  boundary 
layer  equations  are  modified  to  include 
all  pertinent  effects  of  cross  flow  dif- 
fusion and  centrifugal  force.  With  the 
usual  boundary-layer  approximations 

, Re»  1 and  the  reten- 
tion of  cross-diffusion  terms,  required  to 
adequately  describe  boundary  regions  or 
local  shear  flows  formed  near  separation 
plane,  the  Navier-Stokes  equations  can  be 
reduced  to  , 

continuity 

j\  *«*(Pw)i>+f>arx 

x:  momentun4.vov(f?e*X)'^  rj  = 0 
lf)a\i.x-^/3UUy+/>v Uy  + />w  x t±ty 

r,  = - ^47)7V|  r,  a1 


momentum 


(15) 


<f>  momentum  _ 

y>u.  "7  + *r- T 

+/ZUMX  rl  = -M l Vrsm*  P* 

^7+f  h 

energy 

X/>UTt- 


4-  P(  K-l J X [ U,  - jr  -J-  LL  TJ  +-  -LLC  + ^ 

+ -^+t  r^l=0^Wyif'?r7l 

t «/»S*  -J-  j 

umI  r(r-i ) y 4-  ^‘4  ^ 

+ (^H 

The  coordinate  system  is  depicted  in  Pi  . 
(1)  and  the  governing  system  is  net  res- 
tricted tc  small  cross  flow.  Here  the 
boundary  conditions  for  the  boundary  lay- 
er flow  are: 

■>j-0  , 0 } T-T*  Ty  - O 

7-75  , u=ue  , ^ ^ 

III. 2 The  Turbulence  Model 

In  eqs.  (11-15)  it  is  postulated 
that  the  Reynolds  stresses  are  related  to 
the  mean  rate  of  strain  via  a turbulent 
eddy  viscosity;  i.e. 

~ molecular  viscosity 

^tx  2 

A simple  eddy  viscosity  model  is 
based  on  Prandtl's  mixing  length  hypo- 
thesis. For  a three-dimensional  boundary 
layer,  it  is  assumed  that  is  a scaler 
function  independent  of  the  coordinate 
direction.  Accordingly,  the  eddy  viscos- 
ity^) can  be  written  as : 

where  i=(u'+\ V'/1 

Z-  X 8 ta,n/>  ( 44  ~j?~) , * = 0 7 

D = Van  Driest 's  damping  function  = /— 

exfj-y/A* J,  y+= 

A+=  26  tft/piS'/yt'u) 

For  thick  turbulent  boundary  layers 
with  transverse  curvature  effects, 
Cebeci(33)  has  suggested  the  following 
modification  in  the  wall  region 

u'a, 


Recently  four  Independent  experi- 
ments in  different  laboratories  have  found 
that  the  ratio  o*-^  = t /yU  t is  not 

unity  as  assumed  In  the  isotropic  model. 

In  other  words,  the  eddy  viscosity  is  a 
tensor  Instead  of  an  invariant  scalar 
quantity  . For  example,  Bissonnett  and 
Mellor''a)  are  able  to  demonstrate  that 
the  eddy  .iscooity  is  a scalar  only  in 
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the  small  region  near  the  wall,  while  cA* 
decreases  to  an  average  value  of  0.7  in  1 
the  outer  portion  of  the  boundary  layer. 
Three  sets  of  Independent  European  exper- 
imental measurements (35-37)  have  also  sug- 
gested the  value  of  ch-.  can  be  as  low  as 
0.4  through  the  outer  region  of  the  bound- 
ary layer. 

Herein,  the  turbulence  model  will  be 
modified  to  incorporate  the  nonlstrcpic 
form;  the  magnitude  is  dependent  on 

the  coordinate  direction  so  that 

& >^A***U, 

Jt,  = A,  & tan  A JL) 

t*-**  (XT  it)  > 1 

In  this  formulation,  it  has  been  shown 
that  near  the  wall^  , = JL  „ = 0.41  y, 
so  that  c*.  . = 1 is  preserved;  the  two 
eddy  viscosities  will  differ  in  the  outer 
wake  regionj 

In  higher  order  theory,  the  eddy 
viscosity  is  assumed  to  be  proportional  to 
the  transport  properties.  For  example,  in 
Jones  and  Launder ' s model,  is  deter- 

mined by  the  local  values  of  the  density, 
turbulent  kinetic  energy  k,  and  the  dis- 
sipation function,  £ . The  governing  ,-.o\ 

differential  equations  for  K and  € are:^0' 

[U> 

f-*  (“>  * U9) 

WuL)  - Ct />£  Z* 

- C3/>  /tV"',  C,  =/.  +S-,  Ct-z{ i- 0.3 e -/P*1) 

yf»  = ^Al 

Ci  -0.01  etp\-Z.S/[ /+  Rt/S’OjLJ’tf=/j  J1-/.3 

This  model  contains  five  empirical 
constants  which  are  determined  from  ex- 
perimental data.  The  advantage  of  the  k-£ 
model  is  that  the  flow  history  is  taken 
into  account.  With  some  adjustment  on 
the  empirical  constants,  flow  relaminariz- 
ation  or  even  transition  can  be  considered. 

The  following  boundary  conditions 
are  imposed  on  the  equations  for  k and£  : 

y = o,  K=e=  o 

y.S  . u. 

The  turbulent  thermal  conductivity .is 
assumed  to  take  the  following  form;(3oa) 


8nd  * = Cj^JPrt 
Prt=  0.1S--0.4S- 

Finally,  Dhawan  apd  Narasimha's 
intermittency  factor' y t)  is  used  to  mod- 
el the  flow  in  the  transitional  region. 
The  beginning  of  transition  can  be  either 
specified  in  the  (^,  $)  plane  cr  baaed 
on  the  (Re.)  =f  (M  ; criterion.  ( - 

For  the  present  investigation,  we  are  in- 
terested in  evaluating  the  applicability 
of  these  simple  closure  models  for  the 
solution  of  three-dimensional  boundary 
layer  flows,  including  cases  with  cross 
flow  reversal. 


III. 3 Numerical  Methods 

v 

Numerical  computations  in  blunt  ncc- 
region  are  initiated  at  the  stagnation 
point.  A shooting  method  is  used  to 
integrate  eq . (13).  Then  Krause's 
scheme(26)  is  employed  to  solve  eq.  fll'1 
in  a downstream  marching  fashion.  The 
nonlinear  terms  in  the  finite  diff« ren"1 
equations  are  linearized  by  the  Nev.ton- 
Raphson  procedure; 

<"/t  - 2 <ufUa/lr 

where  L denotes  the  Iteration  number. 

The  iteration  continues  until  a 
specified  convergence  criteria  is  at- 
tained. The  solution  is  obtained  with 
an  efficient  "block  tridiagonal"  algor- 
item  described  in  Reference  (40). 

A predictor-corrector  scheme^7’  ' 1 
is  used  to  continue  the  three  dimensional 
boundary  layer  calculation  into  the  super- 
sonic afterbody  region.  Some  modifica- 
tions are  made  in  the  original  P/C  form- 
ulation in  order  to  improve  its  effic- 
iency and  numerical  stability  restric- 
tions. For  turbulent  flows,  it  bar. been 
found  that  a modified  difference'2*!  ) re- 
presentation for  the  lateral  derivative 

, is  superior  to  the  standard 
central  difference,  l.e. 


2A4> 


w >o 


mi 

(21) 


2AQ 

where  M and  K denotes  the  Indices  in  the 
£ and  ^ direction  respectively. 


Kfc  T = -Dy'A'  This  procedure  will  result  in  a ma - 

y / trix  system  which  Is  always  diagonally 

dominant.  Furthermore,  it  eliminates 
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wiggles  which  sometimes  appear  near  the 
cross  flow  separation  line.  At  the  same 
time,  the  numerical  accuracy  of  this  met- 
hod remains  of  second  order. 

In  most  cases,  calculations  were 
initiated  at  the  windward  plane.  With 
the  formulation  given  in  eq.  (2l),  a 
linear  stability  analysis  shows  that  the 
P/C  scheme  is  unconditionally  stable  for 
w > 0,  with  a CFL  condition  on  the  cross 
flow  velocity  w required  if  w<0. 

Because  of  the  various  length  scales 
appearing  in  the  turbulent  flow  (e,g,  the 
laminar  sublayer,  wall  and  wake  region), 
it  is  usually  necessary  to  impose  a co- 
ordinate transformation  or  to  adopt  a 
variable  grid  system.  A simple  mapping 
of  the  following  form  will  serve  this 
purpose : 


y~ln  7 or  y=Y 


When  a variable  mesh  system  is  used,  the 
conventional  three-point  finite  differ- 
ence quotient  can  be  written  as: 


1 1 

Aflll 

M 

\ 

L *' 

( 


|| =*.  CJ)  u3  + CD  U„+  Zc  (rmr 


(a^-a^)+0(a^ 


4 JL 


where 


?4CJj=A^|A^(Ay,  + A^| 

^ ^“A3‘/[A3,(A|+A3t)| 

z*.  */H*(*x<  +A^j 

O)--  i 

Zy  Q)~  2/ [ ajiliiji  FAijjJ 


The  truncation  error  for  U is 
0(  Ay2)  when  uniform  grid  is 
specified,  but  degenerates  to 
0(  a yi  - A y2)  with  a variable  mesh. 
It  should  be  noted  that  U is  always 
second  order  even  for  non-uniform  grids. 

In  order  to  reduce  the  truncation 
error,  one  can  evaluate  U from  the 

governing  equations  and  substitute 

the  resulting  expression  into  the  dif- 
ference quotient,  eq . (22).  In  doing  so, 
a second  order  accurate  system  is  obtain- 
ed for  U.„..  To  illustrate  this  point. 


consider^the  simplified  X 
equation  near  the  wall. 

Liar  expression 


momentum 


(22a) 


One  can  estimate  U by  differen- 
tiating eq.  (22a)  and  substituting 

Into  eq  „ (22),  The  result  is 

u„- &mu, + 2t to u„ ,7, mu, - UAr ~o;  % ^ 

+~oc»f) 

Here  the  term  is  replaced  by  a 

first  order  relation,  J . . 

+ 2^  MM/ 

■/-  0( A y,  - A 

Eq . (23a)  results  In  1 second  rd  1 
accurate  representation  for  U . Most  of 
the  results  presented  here  are^obtalned 
With  this  type  of  formulation. 

An  alternating  direction  Implicit 
methodf^1)  which  Is  unconditi  nally  sta- 
ble for  linear  systems  was  al  ><  mini  l. 
However,  satisfactory  results  wen  ob- 
tained only  for  laminar  flows. 


III.^  Results  fr  r Three -Dime  r.."  ' < na  1 
Viscous  Flow  Calculation:- 

Figures  (6)  and  (7)  depict  the  heat 
transfer  distribution  on  pointed'1*')  and 
blunt  cones. Laminar  and  turbulent 
flows  have  been  studied  and  comparisons 
with  experimental  data  are  encouraging. 
The  flow  properties  at  the  leeward  plane 
of  a sharp  cone  at  large  anglr-  of  att-  :ck 
are  also  predict'  l 1 enably  well 
(Fig.  6) . 

The  Nusselt  number  distribution. on  a 
cone-cylinder-flare  configuration''-’"  is 
shown  in  Fig.  (8).  Significantly,  the 
calculations  do  predict  the  laminar  heat- 
ing o/ershoot  in  a region  whore  a -trong- 
ly  favorable  pi  ure  gradient  exists. 

It  is  noted  that  when  the  flow  undergoes 
a rapid  expansion  (e.g.  at  the  junction 
of  cone  and  cylinder),  the  solution  1: 
very  sensitive  to  the  outer  boundary  con- 
ditions, The  boundary  layer  thickness 
chai  ges  rapidly  and  it  is  no  longer  ad- 
equate to  impose  the  condition  U— ►U 
y— *■&,  An  alternate  prr  si  lure, 
gested  by  Ackerberg  and  Philip-  ) proves 
to  be  satisfactory.  They  use. the  fact 
that  the  velocity  profiles'11’)  exhibit 
the  following  variation*  for  J >5>/ 


u=iU  +47*  exp\-(7-t>)*/ A \ 

U*dt 


(2l>) 


i**  lml 


press  16ns  can  be  written  frr  the  temp'-rature . 
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where  a,  b and  k are  constant  at  f ixed £ 
These  values  are  determined  by  fitting 
eq,  (24)  to  the  velocity  profile  at  the 
outer  edge  of  the  boundary  layer.  This 
procedure  is  incorporated  into  the  im- 
plicit algorithm  described  previously. 

The  numerical  results  for  the  vel- 
ocity and  temperature  profiles  on  a sharp 
cone  are  given  on  Fig.  (9).  The  body 
geometry  and  free  stream  values  cor- 
respond to  Rainbird's  experimental  condi- 
tions. I1*®)  When  the  eddy  viscosity  is 
treated  as  an  invariant  scalar  (i.e  eq. 

(17) ),  agreement  between  the  numerical 
prediction  and  Rainbird's  data  is  good 
for  _s  135°  (Pig.  9)-  But  the  com- 
parison deteriorates  somewhat  as  the  lee- 
ward plane  is  approached.  Similar  ob- 
servations apply  for  the  limiting  stream- 
line inclination  (Fig.  10).  It  is  signi- 
ficant that  the  Prandtl  mixing  length 
theory,  with  a scalar  eddy  viscosity  im- 
plies an  attached  boundary  layer;  the 
experimental  measurements  indicate  that 
cross  flow, separation  occurs  for  $>  l6C° . 
To  the  authors'  knowledge,  this  is  the 
first  time  that  theoretical  results  have 
been  reported  for  the  turbulent  flow  near 
the  leeward  plane  where  secondary  flow 
reversal  has  occurred. 

The  calculated  windward  plane  results 
using  the  two-equation  model  for  Reynolds 
stress  closure  are  also  shown  on  Fig.  9& • 
The  agreement  with  the  data  is  good. 

There  is  no  apparent  advantage  in  using 
this  higher  order  theory,  since  the  sim- 
ple modified  mixing  length  model  leads 
to  an  equally  accurate  prediction  at  the 
windward  plane.  At  the  present  time  num- 
erical computations  using  the  two-equation 
model  have  been  made  only  at  the  sym- 
metry plane . 

In  order  to  improve  the  comparison 
between  the  numerical  results  and  the 
experimental  data  near  the  leeplane, 
some  modifications  of  the  eddy  viscosity 
formulations  are  required.  As  discussed 
previously . experimental  measure- 
ments (3^  -37)  suggest  that  the  eddy  vis- 
cosity is  not  an  isotropic  scalar  quant- 
ity in  the  outer  wake  region.  Although 
the  distribution  ofy*f  and  . 

C1  C2 

across  the  boundary  layer  cannot  be  meas- 
ured accurately  at  the  present  time,  the 
mean  value  of  tne  ratio  ct,  = ^ t2/^tl 

has  been  found  to  be  as  large  as  0.7 
(Ref.  33)  and  as  low  as  0.4  (Ref.  34,  36, 
37).  Fig.  (10)  shows  the  numerical  re- 
sults for  the  surface  limiting  stream- 
line inclination  for  nonisotropic  eddy 
viscosity  distributions  as  given  in  eq. 

(18) .  For  this  calculation  A , = 0.09, 

\ ? = 0.064  and  ol,  = 0.5  are  assumed. 

For  fa > l60°  , there  does  not  appear  to 
be  any  significant  improvement  over  the 
isotropic  model. 

(R3)  (54) 

Bradshaw, Baker  and  Jones w ' 


have  pointed  out  that  the  mixing  length 
theory  may  become  inadequate  to  predict 
the  flows  in  which  the  boundary  layer 
thickness  grows  rapidly  (such  as  flow 
under  adverse  pressure  gradients).  In 
this  case,  because  of  the  convective 
transport  of  turbulence,  the  magnitude  of 
the  mixing  length  in  the  outer  part  of 
the  boundary  layer  does  not  increase  - 
fast  as  the  boundary  layer  thickness.  Ir. 
other  words,  the  eddies  which  may  have 
originated  near  the  windward  plane,  carry 
some  of  the  character  of  the  boundary 
layer  at  an  earlier  stage  of  its  develop- 
ment. Near  the  leeplane  a rapid  thick  n- 
ing  of  the  boundary  layer  occurs . The 
mixing  length  is  in  fact  indicative  of 
the  flow  some  distance  upstream,  where  5 
was  appreciably  thinner.  Consequently, 
the  apparent  value  of  A .in  eqs . (17) 
or  (18)  falls.  Here  an  adjustment  for 
this  effect  is  made  by  reducing  the  con- 
stant X . from  0.09  to  0.063  for  the 
flow  riearthe  leeward  plane  (ip  l45c). 
Results  are  given  on  Figs.  (9E)  and  (1C). 
Crossflow  separation  is  predicted  some- 
what more  accurately  with  this  modifica- 
tion, although  the  improvement  is  only 
marginal . 


IV.  The  Coupled  Viscous  and  Invlscid 
Flow  Computations 

Solutions  for  the  inner  and  outer 
layers  are  now  coupled  in  order  to  take 
into  account  the  visccv s-inviscid  inter- 
action. For  the  lnviscid  (outer)  flcy. .. 
computations,  the  viscous  displacem- ni- 
ls included  by  considering  the  effective 
body  shape  as  modified  by 

reFf(^4>)  = r„(F,4>)*-A  case 

A is  the  three  dimensional  displace- 
ment thickness,  which  must  be  obtained 
from  the  partial  differential  equation: 

^u-s})|=o(26! 

where 

The  body  entropy  value,  " (x),  ir 

estimated  by  the  steamtube  replacement  at 
the  windward  plane,  i.e. 

5„  = 5s(rs  ,o) 

The  procedure  for  evaluating  ? is  dem- 
onstrated in  Fig.  1.  Our  formulation  by- 
passes the  thin  "entropy  or  vortical  lay- 
er" effect  which  can  lead  to  numerical 
instability  when  Z/R^  $>  1. 

For  the  inner  layer  computations, 
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the  boundary  layer  edge  properties  are 
obtained  by  interpolating  the  inviscid 
properties  at  r = r + dcos  0.  Here  d 
is  found  from(4$) 

J=S-&  (27) 


When  the  external  inviscid  flow  is  highly 
rotational,  then  the  following  criter- 
ion(^9a,  5l)  is  used  to  define  the  local 
value  of  S , 


(h  = 0.99 

This  condition  replaces  the  conven- 
tional procedure  of  setting  U = T =0 
as  y— »■  00  . However,  for  adiaSatic^walls 
it  is  necessary  to  revert  to  the  conven- 
tional U =0  condition  to  locate  the 
boundaryylayer  edge . This  implies  that 
the  inner  flow  calculation  contains  part 
of  the  inviscid  vortical  layer. 


The  procedure  for  the  numerical 
computation,  starts  with  the  calculation 
of  the  outer  flow.  The  inviscid  inform- 
ation is  input  into  the  inner  layer  cal- 
culations. However,  when  the  Euler  equa- 
tions are  integrated  at  2—2/  , the 

local  value  of  A is  still  unknown  as 
the  boundary  layer  calculations  have  not 
reached  2 ..  Herein,  a "cyclic  itera- 

tion" procedure  is  employed.  We, initi- 
ate the  calculations  of  the  outer  layer 
by  estimating  two  sets  of  values  A /&  x . 
These  are  obtained  from  Taylor  series 
extrapolations : 


4-  o.S 


o 


and 


am+i  - am-1-  ^ A j A X 


After  completing  the  inviscid  flow 
computations,  this  information  is  input 
to  the  three  dimensional  boundary  layer 
computations  from  which  corrected  values 
of  ||JC  = ( &a/$>x)m  + 1 are  obtained. 

Then  a new  estimate  for^A  /A  x is 
written  as(55>  50  ) 


certain  flow  conditions  and,  therefore, 
becomes  impractical.  However,  a Newtcn- 
Haphson  procedure  (i.e,  eq . (29))  work,, 
satisfactorily  in  all  cases  tested  so 
far.  Usually  3 iterations  are  necessarv. 


The  numerical  procedures  of  coupling 
the  inner  and  outer  flow  have  been  ap- 
plied to  a number  of  test  cases;  howev  r, 
for  the  present  paper  only  one  case  will 
be  presented,  .'’ample  results  consider-:: 
here  are  for  the  supersonic  flow  over  • 
sphere -cone  (0  = 9 ) at  10c  angle  of  at- 
tack. Free-stream . flow  properties  are 
based  on  Widhopf 1 s (^° ) experimental  con- 
ditions. The  inviscid  blunt  body  result 
are  obtained  by  the  time  dependent  met;  - J 
without  the  viscous  displacement  cor- 
rection, since  this  effect  is  small  in 
the  blunt  nose  region  for  the  high  Rey- 
nolds number  flow  (Re0  =1.8  x lC7/ft., 
Mm  = 5)  studied  here.0  The  Blottner  and 
Ellis  formulation  is  used  for  the  inner 
boundary  layer  calculations  which  arr 
initiated  at  the  stagnation  point. 


In  the  supersonic  conical  section 
the  inner  and  outer  flows  are  determined 
by  a marching  procedure.  Viscous  dis- 
placement and  entropy  swallowing  srr  in- 
cluded In  the  Iterative  matching  of  th 
two  regions.  A three-point  Lagrange's 
formula  is  used  to  interpolate  the  bound- 
ary layer  edge  conditions  from  the  in- 
viscid flow  properties. 


The  surface  coordinates  for  the 
blunt  body  boundary  layer  calculations 
are  illustrated  In  Fig.  (11).  The  heat 
transfer  distribution  on  the  body  is  de- 
picted In  Fig.  (12).  Here  the  transiticn 
points  are  specified  from  the  experiment- 
al measurements.  It  Is  shown  in  Fig.  12 
that  the  heating  (for  ? /fo„  . < 5 ) will  be 
underpredicted  by  at  most  ♦'"when  en- 
tropy layer  swallowing  effects  are  not 
Included.  These  effects  will  become 
more  important  at  higher  Mach  numbers 
and  lower  Reynolds  number.  It  is  also 
found  that  the  viscous  displacement  ef- 
fects have  a small  influence  on  the  In- 
viscid pressure  distribution  for  the 
flow  conditions  under  investigations 
(Fig.  13).  Perhaps  this  is  expected 
since  A/r  is  always  less  than  0.1. 
Comparisons  wbetween  the  numerical  re- 
sults and  Widhopf 's  data  (hest  transfer 
and  surface  pressure)  is  good. 


fcA  _ 
b X 


in  (2)  cd 

i vf  - € t 


(2) 


(29) 


The  iteration  procedure  continues  until 
convergence  is  achieved.  During  the 
development  of  the  numerical  program  It 
was  found  that  a fixed-point  iteration 
(i.e.  using  the  most  recently  calculated 
^A/<J  x for  the  next  cycle  of  cal- 
culation) converges  quite  slowly  for 


V.  Summary 

A method  has  been  developed  for 
treating  the  viscous  flow  over  a blunt 
or  sharp  body  at  angle  of  attack.  Here- 
in, a two  layer  model  is  suggested.  The 
inner  region  consists  of  the  three  dim- 
ensional boundary  layer  and  boundary 
region.  Laminar  and  turbulent  flows  are 
considered.  The  governing  system  can 
handle  problems  with  cross  flow  revers- 
al and  Is  integrated  by  predictor- 
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corrector  or  alternating  direction  im- 
plicit method.  The  outer  inviscid  flow 
is  computed  with  MacCormack's  two  level 
finite  difference  scheme,  with  the  bow 
shocks  treated  as  discontinuities.  The 
salient  features  for  matching  these  two 
regions  include  the  effects  of  the  three 
dimensional  viscous  displacement  and 
entropy  layer  swallowing.  It  is  only 
necessary  to  input  the  body  geometry, 
free-stream  flow  properties  and  the  sur- 
face conditions.  The  numerical  results 
include  the  aerodynamic  coefficients, 
heat  transfer  and  the  detailed  flow  pro- 
files . 

The  general  treatment  of  the  prob- 
lem and  the  method  of  solution  are  ver- 
ified by  the  good  agreement  obtained  be- 
tween results  from  the  present  formula- 
tion and  the  experimental  data.  It  is 
observed  in  our  preliminary  results 
that:  (l)  the  simple  scalar  mixing  len- 
gth theory  for  the  Reynolds  stress  ex- 
hibits minor  defects  in  regions  with 
cross  flow  -separation.  Some  adjust- 
ments are  necessary  in  order  to  obtain 
a better  comparison  with  experimental 
data;  (2)  for  numerical  results  not 
shown  here,  the  viscous  displacement 
effects  may  become  more  pronounced  in 
laminar  than  turbulent  flow  and  (3)  the 
entropy-layer  swallowing  is  of  only  min- 
or importance  for  the  examples  consider- 
ed here;  nevertheless,  it  is  expected 
that  this  phenomena  can  become  dominant 
at  hypersonic  speeds  and  for  low  Reynolds 
numbers . 

Work  in  progress  includes  the  fol- 
lowing investigations:  (i)  supersonic 
flow  over  cones  at  large  angles  of  at- 
tack ( d. /©  1.5;  attention  will  focus 

on  the  existence  of  inviscid  solutions), 
(ii)  optimization  of  the  numerical  pro- 
• cedures  for  coupling  the  inner  and  outer 

flow  and  for  the  interpolation  of  the 
boundary  layer  edge  properties,  and 
(ill)  the  application  of  improved  Rey- 
nolds stress  modelling  for  the  three- 
dimensional  boundary  layer. 
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Fig.  6a  Laminar  Heat  Transfer 

Distribution  on  a Pointed 
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Fig.  6b  Turbulent  Heat  Transfer 
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Pig.  10  Surface  Limiting  Streamline 
Inclination 
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Pig.  11  Surface  Coordinate  on 
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Fig.  12a  Heat  Transfer  Distribution 
on  Sphere-Cone  at  10°  Angle 
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Fig.  13a  Surface  Pressure  Distribution 
on  Sphere-Cone  at  10°  Angle 
of  Attack 
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Fig.  13b  Surface  Pressure  Distribution 
on  Sphere -Cone  at  10°  Angle 
of  Attack 
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